Evolution of robust network topologies: Emergence of central backbones 
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We model the robustness against random failure or intentional attack of networks with arbitrary 
large-scale structure. We construct a block-based model which incorporates — in a general fashion 
— both connectivity and interdependence links, as well as arbitrary degree distributions and block 
correlations. By optimizing the percolation properties of this general class of networks, we identify 
a simple core-periphery structure as the topology most robust against random failure. In such 
networks, a distinct and small "core" of nodes with higher degree is responsible for most of the 
connectivity, functioning as a central "backbone" of the system. This centralized topology remains 
the optimal structure when other constraints are imposed, such as a given fraction of interdependence 
links and fixed degree distributions. This distinguishes simple centralized topologies as the most 
likely to emerge, when robustness against failure is the dominant evolutionary force. 



One of the most important characteristics of large net- 
worked systems is their capacity to function reliably. Per- 
haps the simplest paradigm employed in the study of 
robustness of network systems is the theory of percola- 
tion [I], which describes the conditions for the formation 
of a system-spanning connected component. Arguably, if 
a system is not connected to begin with, it is unlikely to 
function properly, regardless of its intended purpose. Re- 
cently [2j, it has been shown that if the interdependence 
of the many elements is taken into account, the random 
failure of a fraction of the system can cause a catastrophic 
outcome, where the relative size of the connected cluster 
falls discontinuously to zero, which otherwise would have 
been a continuous transition, in the absence of interde- 
pendence. 

In this Letter, we investigate analytically and numer- 
ically the most fundamental large-scale structures [3] 
which results in robustness against random failure and 
malicious attack. By constructing a general block-based 
model, and analytically deriving its percolation proper- 
ties, we obtain the topological configuration which opti- 
mizes a well-defined robustness criterion. In the case of 
random failure, we find that a remarkably simple core- 
periphery (CP) topology |U |5] emerges as the most op- 
timal, when no other constraints are imposed other than 
the overall cost in realizing the system. This topology 
remains optimal when one introduces interdependence 
links, and it can entirely suppress the catastrophic failure 
present in fully random topologies. However, in the case 
of targeted attacks, the fully random configuration turns 
out to be the most robust. We also consider the scenario 
where degree constraints are present, and find that the 
bimodal core-periphery structure is replaced by strongly 
dissortative and assort at ive block topologies, for random 
and targeted attacks, respectively. 

Here we consider a stochastic blockmodel [6j[7], which 
defines an ensemble of networks composed of B discrete 
node blocks, where n r is the number of nodes in block 
r and e rs is the number of edges between blocks r and 
s (or twice that number if r = s). Each block is al- 
lowed its own degree distribution, p r k . We also consider 



interdependence edges by defining a distinct matrix e rs 
and degree distributions p r y in an analogous fashion. As 
in [2 J, if a node u is dependent of another node v, u fails 
automatically when v fails, and vice-versa. In the case of 
multiple dependencies [8j, all support nodes must fail in 
order for the dependent node to fail (i.e. they are redun- 
dant). Assuming that the values of n r are large enough, 
this ensemble becomes locally tree-like [9 J , thus it is pos- 
sible to adapt the epidemics-based generating function 
formalism [10-12J, which becomes exact in the limit of 
large networks. If we define Uf as the probability that a 
node belonging to block r is not in a macroscopic com- 
ponent via one of its neighbors, the dilution variable <j) r k 
as the fraction of nodes belonging to block r and with 
degree k which are not removed from the network, and 
additionally the interdependence dilution <p r denned as 
the fraction of nodes from block r which were not re- 
moved due to the failure of the support nodes, we can 
write the following self-consistency equations (see Sup- 
plemental Material for a derivation), 
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where m rs = e rs /n r n r is the asymmetric matrix defin- 
ing the fraction of edges adjacent to vertices of block r 
which are also adjacent to block s, where n r = e rs /n r 
is the average degree of block r. S r = (j) r S® is the frac- 
tion of nodes of block r which belong to a macroscopic 
component, and S® = /o (1) — fo( u r)i given as a func- 
tion of the diluted degree generating function f${z) = 
EfcPfc^fe**- Furthermore, f{{z) = £ fe qk^l +1 z k gen- 
erates the diluted excess degree distribution of block r, 
with q k = p r k+1 (k + 1)/ n r . Note that we have f{(z) = 
f r ' (z)/g r ' (l) where = ^2 k p 7 k z k generates the undi- 
luted degree distribution of block r. The generating func- 
tion /o(z) = ^2%p r k z k describes the degrees correspond- 
ing to interdependent edges alone. The total fraction of 
nodes which belong to a macroscopic component is sim- 
ply S = s }2 r w r S ri where w r = n r /N is the fraction of 
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FIG. 1. (a) Block structure with 10 blocks of equal size. Each 
node represents a block, and black (red) edges correspond to 
connectivity (interdependence) edges. Squares (circles) have 
exponentially (Poisson) distributed connectivity degrees with 
Ki — 2.25 (Ki = 1.75), and ki = 8 (ki = 6). All interdepen- 
dence degrees are exponentially distributed, (b) The matrices 
e rs (black squares) and e rs (red squares) that correspond to 
this ensemble, (c) S a function of <fi, for random failures. The 
solid lines are results obtained via Eqs.[l]and[2] and the filled 
symbols are obtained empirically with a network realization 
of size N = 10 6 . The two discontinuous jumps, correspond to 
the outer and inner rings failing in sequence. 

nodes belonging to block r, and the total dilution in the 
network is likewise (j) = Y^ r ,k w rP r k^k P3]- 

Although quite straightforward, this parametrization 
is a generalization of many scenarios considered in 
the literature, such as two fully interdependent net- 
works [2j, the "network of networks" topology [T4UT6] . 
single networks with mixed connectivity and interdepen- 
dence edges [T7HT9] , as well as networks with only con- 
nectivity edges and arbitrary degree distribution and cor- 
relations [20J, all of which correspond to specific choices 
of the parameters defined here. In Fig. [T] is shown an 
example of a block structure with both support and in- 
terdependence links, as described in the legend. It also 
shows a comparison with the empirical percolation pro- 
file of networks realized from this ensemble, which is in 
excellent agreement with the theoretical prediction. 

We are interested in obtaining the block topologies for 
which the robustness against failure or attack are opti- 
mal. Following [21 J we will consider the total robustness 
of a network ensemble to be given by 

R = 2 f 5(0)#, (3) 
Jo 

where the factor 2 is chosen so that R G [0,1], with 
R = 1 being possible only in the hypothetical scenario 
S(4>) = (/>, attainable only for infinitely dense networks. 
Thus our ultimate task is to find the parametrization of 
the blockmodel ensemble which maximizes Eq. [3J under 
suitable constraints. In the following, we focus on the 
special case with uniform failure of nodes within a sin- 
gle block, (jf k = <j) r \ hence /q(z) = (j)r9o( z )- The total 
dilution becomes simply <\> = J2 r w r <p r . [Without loss of 
generality, since heterogeneous degree-based dilution can 
still be achieved if nodes of different degrees belong im- 
plicitly to different blocks.] In the case of random failure 
we have simply, <j) r = <j). In the case of targeted attacks 
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FIG. 2. (Color online) (a) S as a function of 0, for random 
failures, different (3 values, (k) = 2 and (k) — 0. The solid 
lines are results obtained via Eqs. [I] and [2] and the filled 
symbols are obtained empirically with a network realization of 
size N — 10 6 . The dashed curve close to (j) — 1 corresponds to 
targeted attacks for (3 — » oo; (b) Correlation matrix e rs /e r e s 
(top), block sizes w r and average degrees K r , for different (3 
values; (c) Ensemble samples with (k) = 2 and N = 10 3 nodes 
for different (3 values. The orange (square) nodes belong to 
the core block. 

we will use (j) r oc e K r(i-b)/b^ w \ iere 5 ^ [q, 1] must be so 
chosen to achieve a desired total (j). 

Instead of directly finding the network topology 
parametrization which maximizes R, we proceed by ob- 
taining the configuration which delivers a specified value 
of R, and is otherwise maximally random. We do so in 
order to identify network structures with a varied de- 
gree of robustness, and to isolate the most fundamen- 
tal characteristic responsible for its increase. Hence, we 
consider a null model of robust networks, where superflu- 
ous characteristics are explicitly discarded. More specifi- 
cally, we maximize the entropy of network ensemble [22J , 
subject to the constraint that the average robustness is 
fixed [23J. This amounts to finding the critical points of 
the Lagrangian A = E — j3(R — R*), where E = lnfi 
is the microcanonical entropy (with Q being the number 
of different network realizations for a specific ensemble 
parametrization), R and R* are the actual and desired 
average robustness respectively, and /3 is a Lagrange mul- 
tiplier. Instead of fixing R*, one may fix j3, and this 
becomes equivalent to minimizing the free energy of the 
ensemble, 

T = -NR - E//3, (4) 

since A = ^(J^/N + R*), and f3 = —fi/N. In Eq. [I| the 
value — NR plays the role of average energy, and f3 is 
the inverse temperature, which can also be interpreted 
as a "selective pressure" [24J, since for j3 = the en- 
tropy is strictly maximized (i.e. all networks occur with 
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equal probability), and conversely for (3 — ^ oo the av- 
erage robustness is strictly maximized. In our case we 
have E = S + <S, where S is the stochastic blockmodel 
entropy [25J, 



Pk( ln Pk + lnfc!) - ^^Z e ^ ln ( ~T" 



(5) 

which includes the entropy of the degree distributions of 
the individual blocks, and S is equivalent ly defined as a 
function of e rs and p 7 ^ [26J . 

Solving the system given by Eqs. [I] and [2] cannot be 
performed analytically, and thus the same holds for ob- 
taining R and T . Hence we must resort to solving Eqs. [I] 
and [2] numerically (by simple iteration) , in order to ob- 
tain Eqs. [3] and [4] The minimization of Eq. [4] can then 
be performed within arbitrary precision with any suit- 
able minimization algorithm (see Supplemental Material 
for details). 

The ensemble can have additional constraints, which 
must be imposed when minimizing T . One such con- 
straint we will consider throughout this paper is that the 
average degree (k) will be fixed at a given value, which 
represents the putative cost of adding extra edges to the 
network, compared to simply replacing them. Further- 
more, in order to restrict the number of degrees of free- 
dom of the model to a manageable size, we will assume 
that all blocks have the same type of degree distribution, 
which can be different only in their average value, n r . 
This does not alter very strongly the type of networks 
which are ultimately attainable, since in principle one 
can construct many different topologies by composing a 
larger number of such blocks. We will also forbid nodes 
with degree zero, since these never belong to a macro- 
scopic component. The degree distribution we use is a 
modified Poisson distribution p r k = (1 — Sk,o)K>^/(e K " r — 
l)fc! [identical to a regular Poisson, but with support 
k > 1], which is generated by go(z) = ze^- i)(2-i) ? 
which of course implies that n r > 1 for all r. The 
same distribution will also be used for the interdepen- 
dence degrees, Note that this simplifies the entropy 
to S = ^2 r n r ln(e Kr - 1) - \ Y, rs e rs ln(e rs /n r n s ), and 
equivalently for S. Finally, the total number of blocks B 
will also remain fixed during the minimization of T . 

In absence of any further constraints, we start with the 
case of random failures, and no interdependence links 
(k r = 0). We observe that increasing j3 leads from a 
fully random graph (/? = 0) to a topology which can be 
significantly more robust (see Fig. [2^l) . For B = 2 blocks, 
this topology is a core-periphery (CP) structure [4j |5], 
where one block is much smaller and has a much larger 
average degree, K r , and attracts a large portion of the 
edges. One might expect this to be a special case when 
one imposes 5 = 2, and that different robust topologies 
may arise for B > 2. However, this turns out not to 
be the case, and we find exactly the same CP topology 



for any value of B larger than two, in the sense that 
two or more blocks may be merged together, with the 
ensemble remaining equivalent, until only B = 2 blocks 
are left. In Fig. ^2jp are shown two resulting ensembles 
for different j3 and with B = 3, for (k) — 2. The size 
and average degree of the core block become smaller and 
larger, respectively, for larger values of j3. Examples of 
how such networks can look like are shown in Fig. |2]2. 
For sufficiently large /?, the vast majority of periphery 
nodes are connected exclusively to core nodes, which are 
much smaller in number (infmitesimaHy so for j3 — » 00), 
and which are connected among themselves, forming a 
central backbone [27J. 

Observing how the value of R changes with f3 exposes 
a structural phase transition at a critical value of j3 = (3* , 
below which the networks are fully random (e rs oc e r e s 
and K r = (fc)), and above which the CP topology is ob- 
served (see Fig [3^l) . The value of /3* becomes smaller for 
(k) close to two, and larger for either sparser or denser 
graphs. The CP topology is particularly advantageous 
for sparse graphs, with lower (fc), but is less so for denser 
graphs, as can be seen in Fig. [3^i, which shows a compar- 
ison between the f) = (random) and j3 — 00 (optimal) 
cases, as a function of (k). 

When node interdependence is enforced {k r > 1 [28]), 
the same CP structure emerges for larger values of /?, 
as shown in Fig. [4b. For f3 —> 00, the percolation tran- 




FIG. 3. (Color online) (a) Robustness R (Eq.[3| as a function 
of /3 for different average degrees (k), with (k) = 0; (b, c) 
Robustness difference between optimal ensembles (f3 — »• 00) 
and random networks (f3 = 0) with the same (k), for (b) 
(k) = and (c) (k) > 1. 
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FIG. 4. (Color online) (a) Same as Fig.jiji, with (k) = 4, (k) = 
1.1.; (b) Same as Fig. [2J3, and the same constraints as in (a). 
The correlation matrix e rs /e r e s and average interdependence 
degrees k r are shown in red; (c) Same as Fig.[2ji, with (k) = 2, 
(k) = 1.1 and N = 10 4 . The red edges are interdependence 
links. 
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FIG. 5. (Color online) (a, b) Same as Fig. |2^l for (a) random 
failures and (b) targeted attacks. In both cases the values w r 
and 

Kir are kept fixed, (c, d) Ensemble samples with N = 
10 4 nodes, f3 —> oo, for (c) random failures and (d) targeted 
attacks. 



sition changes from abrupt to continuous (see Fig. |4p). 
The average interdependence degrees, k r , become corre- 
lated with the connectivity degrees, K r , with most inter- 
dependence edges lying between two core nodes, which 
therefore benefit more from the extra redundancy, to the 
indirect benefit of the periphery nodes as well. In com- 
parison to the (k) = case, denser networks benefit more 
from the CP topology, since it removes the abrupt phase 
transition. The benefit diminishes for larger (fc), since 
the larger redundancy already provides an inherent ro- 
bustness to random networks. A special case is (k) = 1, 
which shows a lower improvement, because it strictly for- 
bids intra-core redundancy. 

As seen in Figs. |2^i and [4^, the CP topology is ex- 
cellent against random failure, but is extremely fragile 
against targeted attacks, in which the core nodes are re- 
moved before the rest (the value of S vanishes for an 
infinitesimal value of 1 — <fi). Instead, if one optimizes 
R against targeted attacks, one obtains that the random 
block topology is the most robust. The only exceptions 
are very sparse networks with (k) < 2 (see Fig.^). Since 
these networks are close or below the percolation thresh- 
old even for <j> = 1, the CP topology is still an improve- 
ment over the fully random case, even if gets destroyed 
very quickly by a targeted attack. 

In some realistic situations there are degree constraints 
restricting the accessible topologies [21J. Here we intro- 
duce this by forcing that the blocks have prescribed sizes 
and average degrees (see Supplemental material for de- 
tails). As can be seen in Fig. [5^i, in the case of random 
failure, the optimal topology approaches the CP config- 
uration as much as the imposed constraints allow, and 
the inter-block connectivity becomes strongly dissorta- 
tive, with the nodes with largest degrees serving as inter- 



mediaries between the nodes with the lowest degrees and 
an innermost "core" composed of more of such alternat- 
ing layers. For targeted attacks the resulting topology 
is quite different: The blocks with higher hv r are ex- 
pelled" from the network, since they are the first to be 
removed in an attack. The resulting topology is highly 
assortative, with blocks with similar K r connecting pref- 
erentially to themselves. This is qualitatively equivalent 
to the "onionlike" topology found empirically in [5T]. An 
exception to this are the blocks with small n r . These do 
form a CP structure, since an assortative connectivity 
would have very weak percolation properties even when 
no nodes are removed, i.e. = 1 (similarly to the case 
(k) < 2 discussed previously). Very similar results are 
obtained when interdependence is introduced (see Sup- 
plemental material). 

The topologies described above arise when no con- 
straint other than those considered is being imposed, 
and the ensemble entropy is maximized. This is not of- 
ten the case for real systems, since they may be subject 
to competing selective pressures, with robustness being 
only one of them, as well as other topological restrictions, 
frozen accidents, etc. Nevertheless, the analysis pre- 
sented here suggests that a simple core-periphery struc- 
ture is the most natural configuration which achieves 
robustness against random failure (and that a random, 
non-centralized topology works best against targeted at- 
tacks). Therefore, it not a surprise that a CP topol- 
ogy is indeed found to some extent in many real sys- 
tems [4], such as the Internet [29J, and gene regulation 
networks [24J. Also, our results indicate that other fea- 
tures found in real systems, such as scale- free degree dis- 
tributions, are not a strict necessity for robustness, and 
do not arise naturally for this purpose, and thus may 
be simply a byproduct of a non-equilibrium organization 
process. 
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GENERATING FUNCTION PERCOLATION 
FORMALISM 

Following [1-3], we derive self-consistency equations 
for the percolation properties of stochastic blockmodel 
ensembles as follows. If we let u r describe the probabil- 
ity that a node belonging to block r is not in a macro- 
scopic component via one of its neighbors, the dilution 
variable (f) r k as the fraction of nodes belonging to block r 
and with degree k which are not removed from the net- 
work, and additionally the interdependence dilution <p r 
defined as the fraction of nodes from block r which were 
not removed due to the failure of the support nodes, we 
can write, 



u r = ^ m 



+ ■ 



<)ql (1) 



where m rs = e rs / n r hvf is the asymmetric matrix defining 
the fraction of edges adjacent to vertices of block r which 
are also adjacent to block s, and q k = p r k+1 {k + 1) / n r 
is the excess degree distribution of block r, and k t = 
e rs /n r is the average degree of block r. This can be 
rewritten in compact form as 
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where f[(z) = Y^k Qk ( t ) k-\-i z is the- generating func- 
tion of the excess degree distribution. Note that we 
have f[(z) = f r ' (z)/g r ' (l) where g r = T, k Pk zk and 
fo( z ) ~ ^ZkPk^k^ are the generating functions of the 
original and diluted degree distributions, respectively. 
The fraction of nodes of block r which belong to a macro- 
scopic component is given by S r = 4>rS®, where S® = 
/o (1) — fo( u r)- The total fraction of nodes which belong 
to a macroscopic component is simply S = ^ r w r S ri 
where w r = n r /N is the fraction of nodes belonging to 
block r, and the total dilution in the network is likewise 

<t> = T, r ,k w rPk ( t>k- 

In the case of no interdependence edges, we have sim- 
ply cj) s = 1, and Eq. 2 describes S r completely. Oth- 
erwise, an additional fraction of nodes will be removed 
from the macroscopic components, if the nodes to which 
they depend do not belong to the component. The frac- 
tion (j) r of nodes which are not removed in this way can 
be obtained by 
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k 



Pfe[l-(l-E 8 ^-5?)* 
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where k are the interdependence degrees, and p 7 ^ is the 
interdependence degree distribution of block r. Here it is 
assumed that the interdependence is undirected (if node 
u is dependent on v, then v is also dependent on u), and 
that at least one dependency must not fail in order for 
the dependent node not to fail (i.e. redundant interde- 
pendence), and furthermore that a node with k = is 
autonomous (see below for variations of these rules). We 
may write this equation in a more compact form as, 

4>r = 1 " / r (l " E>rA°) + / O r (0), (4) 

where /q(z) = ^2 k p r ^z k is the interdependence degree 
generating function. Note that the right hand side of 
Eq. 4 depends only on u r variables (via 5^), and thus 
may be included directly into Eq. 2, resulting in a to- 
tal of only B equations. Similarly to [3], here we omit 
the detailed description of the failure cascades as done 
in e.g. [2], however these correspond simply to repeated 
iterations of Eqs. 2 and 4, where one must first obtain the 
convergence of the u r variables, and then of the (j) r vari- 
ables, in succession. However, if one is only interested in 
the steady state, simultaneously iterating Eqs. 2 and 4 is 
more straightforward. 

One may also consider different interdependence rules 
by slightly modifying Eq. 4. For instance, in the case of 
non-redundant interdependence, where all dependencies 
are essential, so that if any one of them fails, the de- 
pendent node fails as well, can be described by replacing 
Eq. 4 by 

^ = /o r (E>rA°)+/o r (0). (5) 

Furthermore, one may also consider the situation where 
the interdependence is directed, i.e. if node u depends 
on node v, v does not necessarily depend on u (in this 
situation it assumed that e rs is also directed, and thus 
asymmetric, and the variables k describe the dependence 
in- degree). This can be implemented by an almost trivial 
modification, simply by replacing S® — > S s in Eq. 4 or 5, 
where we have, as before, S r = 4> r S®- This is simply 
because a node will fail either if it does not belong to a 
macroscopic component, or if its dependence neighbors 
fail. However, in the case of directed interdependence, 
the likelihood that two nodes are mutually interdepen- 
dent is vanishingly small, and thus the support node will 
depend exclusively on other nodes, which will themselves 
not fail with probability given by <j) r . This introduces a 
direct coupling between the variables <p r (which is only 
indirect in the undirected case). Although seemingly in- 
nocuous, this modification can drastically change the per- 
colation properties of the system, in particular when the 
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FIG. 1. Robustness R to uniform node removal of random 
networks with one block B — 1 , as a function of (k), and dif- 
ferent values of (k), for (a) undirected interdependence and 
(b) directed dependence; (c) Relative size of the largest com- 
ponent S as a function of the dilution 0, for a random graph 
with (k) = 5, directed dependencies and different (k) values. 
The solid lines are results obtained via Eqs. 2 and 4 (with 
Ss — > S s )j and the filled symbols are obtained empirically 
with a network realization with N = 10 6 nodes. 



average degrees k r are close to one. Indeed if all nodes 
have k = 1, the failure of an infinitesimal fraction of 
nodes will always cause a global cascade, independently 
of how well connected the network is. This is easy to 
understand, since the removal of a single node will neces- 
sarily cause the failure of another node, which will trig- 
ger the failure of a third one, and so on, until a loop is 
reached. The average size of these cascades diverges, and 
thus the values of S r will drop to zero. The robustness of 
random networks (B=l) with different values of (k) and 
(k) are shown in Fig. 1, for both connectivity and inter- 
dependence degree distributions given in the main text 
if (k) >loip k = (k)S kA + (1 - (k))S ki0 if (k) < 1. For 
undirected interdependence, we have that R approaches 
one for (k) — » oo, for all values of (k). However for di- 
rected dependence, we have R dropping to zero as (k) 
approaches one, independently of (k). 

In all cases, the percolation threshold can be obtained 
by examining the Jacobian J = {d(u f r \^ r )/d(u 3 \^) 3 )}, 
where u' r and <\>' r correspond to the right-hand side of 
Eqs. 2 and 4, respectively. If the largest eigenvalue A 
of J at the fixed point u r = 1 and <p r = /o(0) exceeds 
one, this solution is unstable, and thus some value of S r 
must be finite, and hence the condition A = 1 marks the 
percolation transition. 



MINIMIZATION OF THE FREE ENERGY, T 



As mentioned in the main text, the computation of R 
and thus T must be done numerically, by iterating Eqs. 2 
and 4 until convergence within some desired precision is 
reached. This means that the minimization of T must 
also be performed numerically. We used different algo- 
rithms, depending on the set of constraints used. In the 
case of (k) = 0, the L-BFGS-B quasi-Newton method [4], 
with gradient obtained with finite differences, delivered 
fast and reliable results if the number of degrees of free- 
doms were not very large (B — 2 or 3). For larger number 
of blocks, the non-gradient COBYLA algorithm [5] per- 
formed better. In the case (k) > 0, the discontinuity in 
S(<p) makes the gradient estimation unreliable, and the 
COBYLA algorithm is a better choice. In both cases, 
when improved precision was necessary, a repeated se- 
quential one-dimensional minimization over all variables 
using Brent's method [6] was performed, at the expense of 
longer computational time. In the case of very large val- 
ues of an initial estimation using the ALPSO stochas- 
tic swarm-based algorithm [7] was used, which was then 
improved upon, whenever possible, with COBYLA. 

In all cases, the constrained optimization problem 
was transformed into an unconstrained one, via appro- 
priate transformations. The normalization constraints, 
w r = 1 and ^2 rs e rs = N(k) can be imposed simply 
by transforming w r = Cexp(iD r ), and e rs = Z}exp(e rs ) 
with C" 1 = Yj r w r an d D ~ X = ^2rs e rs/N(k), and 
solving for w r and e rs . Note that we can set, e.g. 
w = and e 00 = so that only 2B + B(B - l)/2 - 2 
free variables remain. The degree constraints K r > 1 
can be imposed via the transformations n r = k r if 
ftmin = min({ft r }) > 1 otherwise + 1 where 

a r — (hz r — 1)/^2 S w s (k s — ft m in). The new row and col- 
umn sums of e rs given by e rs = e sr = Nw r n r 
can then be enforced in the e rs matrix via a Sinkhorn 
scaling [8]. The same transformations can be used for 
the e rs matrix. With all these transforms in place, the 
optimization problem becomes fully unconstrained, with 
w r , e rs and e rs lying in the interval [—00,00]. Note that 
the variable N is only included in the text to slightly 
simplify the notation, but it is entirely unnecessary dur- 
ing the optimization, since one needs to work only with 
normalized variables. 

In all minimizations performed, the same solutions 
were found, regardless of the algorithm used or the initial 
configuration. Although the different methods have dif- 
ferent precision and convergence speed, the results were 
always mutually consistent. Additionally, the indepen- 
dence of the results on the initial configuration implies 
the results obtained were always a global minimum of the 
free energy. 
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FIG. 2. (a, c) Fraction S of nodes in the macroscopic com- 
ponent as function of the dilution 0, for different /3 values, 
fixed values of w r and ft r , and (k) = 1.2, for random failures 
(a) and targeted attacks (c). The solid lines are results ob- 
tained via Eqs. 2 and 4, and the filled symbols are obtained 
empirically with a network realization of size N = 10 6 ; (b, 
d) Correlation matrices e rs /e r e s (black) and e rs /e r e s (red), 
block size w r and average degrees n r (black) and k r (red), for 
/3 —> oo, for random failures (b) and targeted attacks (d). 

IMPOSING DEGREE CONSTRAINTS 

As mentioned in the main text, in some realistic situ- 
ations there are degree constraints restricting the acces- 
sible topologies [9]. Here we introduce it by forcing that 
the blocks have prescribed sizes and average degrees, fol- 
lowing w r oc ft~ 7 and n r = e cr , with re [0,5 — 1] and 
c = In ft max / (B — 1). This allows us to have a very broad 
range of degrees, similar to what is found in scale-free 
networks. In Fig. 5 of the main text are shown the re- 
sults for B = 6, 7 = 2 and ft m ax = 10 3 . Very similar 
results are found by varying £?, 7 and K max (not shown). 



Here we expand briefly on these results by introducing 
interdependence, i.e. (k) > 1. We do so in slightly dif- 
ferent way than before, by imposing a separation of the 
network into two parts of equal size, such that connec- 
tivity (interdependence) edges are only allowed between 
nodes in the same (different) part. This is to mimic, e.g. 
the situation of computer and electricity infrastructure, 
where interdependence is only meaningful between the 
two different types of nodes. With this constraint, we 
obtain a configuration very similar to the case without 
interdependence (see Fig. 2), with dissortative and assor- 
tative connectivity for the random and targeted removal, 
respectively. However, the interdependence connectiv- 
ity is slightly more elaborate: For random failures, the 
values of Kj T are correlated with n ri and the e rs matrix 
displays a more assortative connectivity, with blocks de- 
pending on other blocks with similar average degree. In 
the case of targeted attacks, the interdependence shows 
a CP structure for blocks with low degree, but otherwise 
follows the same segregated pattern as the connectivity 
edges. 
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